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ABSTRACT 

We give an explicit formula for the membrane potential of cells in terms of the intracellular 
and extracellular ionic concentrations, and derive equations for the ionic currents that flow through 
channels, exchangers and electrogenic pumps based on simple energy considerations and conservation 
laws. We demonstrate that the work done by the pump is equal to the potential energy of the cell 
plus the energy loss due to the downhill ionic fluxes through the channels and the exchanger. 
Our equations predict osmotic pressure variations. The theory is illustrated in a simple model 
of spontaneously active cells in the cardiac pacemaker. The simulated action potential and the 
five currents in the model are in excellent agreement with experiments. The model predicts the 
experimental observed intracellular ionic concentration of potassium, calcium and sodium. We do 
not see any drift of the values for the concentrations in a long time simulation, instead we can obtain 
the same asymptotic values starting with equal intracellular and extracellular ionic concentrations. 

I. INTRODUCTION 

Ionic models of the cellular membrane potential are typically guided by extensive data sets from voltage clamp 
experiments. The purpose of this paper is not to introduce a new ionic model incorporating the many experimentally 
observed currents (Boyett, 1996), or to critically review other models, but merely to give an alternative formulation 
based on simple principles from classical physics. Instead of trying to reproduce observed voltage clamp data and 
tail currents of activation, deactivation and inactivation, we demonstrate that simple energy considerations and 
conservation laws can be useful when formulating a model. 

We derive equations for ionic currents flowing through channels, exchangers and electrogenic pumps. These are 
based on the Boltzmann distribution law (Boltzmann, 1868), which states that a particle in thermal equilibrium spends 
less time in states of higher energy than in states of lower energy, the Markov assumption (Markov, 1906) which says 
that the transition probabilities to the next state of a stochastic system is only dependent on its present state, and 
the principle of detailed balance (Onsager, 1931) which says that the microscopic laws of physics are invariant with 
respect to the reversal of time. Our equations were inspired by Ehrenstcin and Lecar's model of channel gating (1977), 
MuUins' model of the Na^, Cs?^ exchanger (1977), and Chapman's model of the Na^, K+ pump (1978). 

Since the intracellular ionic concentrations are treated as dynamical variables, we show that the standard differential 
equation for the membrane potential can be replaced with a new algebraic equation. We demonstrate that the 
pump generates a transmembrane voltage gradient, a transmembrane pressure gradient, transmembrane concentration 
gradients, and downhill ionic fluxes through the exchanger and the channels. 

The theory is illustrated with a simple model of spontaneously active cells in the rabbit sinoatrial node. The 
observable parameters in the model are based on the experiments of Shibasaki (1987), Hagiwara et al. (1988), 
Muramatsu et al. (1996) and Sakai et al. (1996). The non-observable parameters in the model are determined 
numerically, in the same way as in an earlier study (Endresen, 1997a), by comparing the action potentials generated 
by the model with the action potentials recorded by Baruscotti et al. (1996). 



II. DERIVATION OF THE EQUATIONS 
A. Nernst Equation 

There are two basic principles behind the average motion of particles. The first is diffusion which is general; the 
second applies only to charged particles such as ions in solutions. Simple diffusion is described by the empirical law 
of Pick: 
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$= -ukTV[S] , 



(1) 



where (p is the ionic flux, [S] the concentration of ions and u the ratio of the velocity to the force acting on a particle, 
known as the mobility. The empirical law of Ohm describes the net motion of charged particles in an electric field: 

= -zeu[S]Vw, (2) 

where z is the valence, e the elementary charge and v the potential. Since we assume that the temperature is constant, 
we can neglect the thermal flux given by Fourier's empirical law. The fact that the mobility in Pick's law must be 
identical to the mobility in Ohm's law was first noticed by Einstein (1905). If we combine equations (|l|) and (|^), the 
total flux of ions due to diffusion and electric forces becomes: 



isi=p(i|;)] . (3) 



The equilibrium potential where the flux is zero can be found by integrating equation from the intracellular (i) to 
the extracellular (e) side of the membrane: 

kT , [S]c 

vs = Vi - I'o = — In -— , (4) 
ze [S\i 

where [S]i and [S]e refer to the intracellular and extracellular concentrations. This equation, first stated by Nernst 
(1888) is based only on the empirical laws of Ohm and Fick and the relation of Einstein. 

The same formula can be derived in a more general way using the Boltzmann factor (Boltzmann, 1868). The 
relative probability at equilibrium that an ion is at the intracellular or extracellular side of a cell membrane is 

Pe [S]c / Aw\ 

- exp , (5) 



p, [S]i ^^^^'V kTj 

where Aw is the energy difference between these intra- and extracellular states , T is the absolute temperature 
and k is Boltzmann's constant. If we consider ions of valence z the energy difference between the intracellular and 
extracellular side can be written: 

Aw = ze{vc — Vi) . (6) 

This together with equation (^ yields equation (^. The equilibrium potentials for the predominant cellular cations 
are then: 

kT. [K]e ,_, 
«K = — In -7—- , 7 
e [K]i 

kT Ca c 
kT [Na]e 



B. Ionic Channels 



We assume that ionic channels are either completely open or completely closed and randomly fluctuate between 
these states. Each channel is here assumed to be a simple Markov process (Markov, 1906), described by first order 
kinetics:: 

^ f O , (10) 

where the rate constants a and /? are functions of transmembrane voltage and control the transitions between the 
closed (C) and the open (O) states of the channel. The rate for a closed channel to open is a, and f3 is the rate for an 
open channel to close. Let x denote the average fraction of channels that are open, or, equivalently, the probability 
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that a given channel will be open. We may say that the ionic flux through an ensemble of channels is regulated by a 
sliding door whose position is x. This yields: 



— = a{l — x) — px = 

at T 



where 
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(11) 
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Here x^o denotes the steady state fraction of open channels and r the relaxation time. Let us assume that the energy 
difference between the closed and open positions is given by 



Aw = Wopcn - Wcloscd = ^(Vx " v) 



(14) 



where q is an gating charge, usually q « 4e, the term qvy^ is due to the difference in mechanical conformational energy 
between the two states and qv represents the change in electrical potential energy due to the redistribution of charge 
during the transition. At equilibrium, dx/dt — 0, and the ratio of the probabilities for a single channel to be in the 
open state or the closed state is, 



\- Xo 



a 



(15) 



This relation is known as the principle of detailed balance (Onsager, 1931). The same ratio is given by the Boltzmann 
distribution (Boltzmann, 1868), 



1 - Xoo \ kT ) 



(16) 



Thus, from equations (|l|), (|l|) and (|T|): 



1 + exp 
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and 
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exp 


(Vx — v) 
kT 
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Therefore, making the ad hoc assumption that the product of forward and backward reaction rates is constant, the 
symmetric choice for a and /3 is: 



a = A exp 
/? = A exp 



(t'x - V) 



kT 
2e 



(I'x - v) 



kT 
2e 



(19) 
(20) 



where A is a constant. Taking A to be constant gives the maximum relaxation time at the voltage where Xoo — 1/2. 
The relaxation time is then: 



1 

a + /3 



= 2A cosh 



qjvi - y) 

2kT 



(21) 



3 



1. Potassium channels 



If the flux of ions is driven by Ohm's law and regulated by the fraction of open channels a;, the membrane current 
through potassium channels becomes: 



iK = gKx{v ~ vk) 

dx _i , I w — 
— — cosh 
dt ^ 



kT 
2e 



1 + tanh 
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(22) 



where tk — 1/2 A is the maximum value of the relaxation time, gy^ is the conductance, is given by equation (|^), 
where we have used the mathematical identity: 



l[l + tanh(0)]= ^^^^^^_^^^ 



(23) 



2. Calcium and sodium channels 



The calcium and sodium channels also have an inactivation mechanism in addition to the above activation mechanism. 
We can view these mechanisms as two independent Markov processes, or equivalently two independent sliding doors, 
which are both affected by voltage. An ion can only go through if both sliding doors are at least slightly open. The 
activation mechanism is very fast, so we use the steady state fraction of open channels (equation (p!8|)) for this. The 
maximum time constant of inactivation for calcium and sodium channels are of the same order of magnitude as the 
maximum time constant of the activation of the potassium channel, thus: 



dt 



.9Ca f{v- WCa) M 1 + ^^"^^ ■ 



^Ca cosh 
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kT 
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and, 
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(25) 



where ijca and t/Na are the conductances of the calcium and sodium currents respectively, vca and w Na are given by 
equations (0) and (m, Vd and are the half-activation potentials, and V{ and Wh are the half-inactivation potentials. 



C. Na+,K+ Pump 

The Na,K-ATPase is found in the plasma membrane of virtually all animal cells and is responsible for active 
transport of sodium and potassium. Low sodium concentration and high potassium concentration in the cytosol are 
essential for basic cellular functions such as excitability, secondary active transport, and volume regulation. In our 
model, the Na+,K+ pump is the only energy source. We shall assume that the following equation is a complete 
macroscopic description of the pump reaction (Chapman, 1978): 

ATP + 3Na+ + 2K+ ^ ADP + Pjo + 3Na+ + 2K+ , (26) 

where ATP, ADP and Pjo are adenosine triphosphate, adenosine diphosophate and inorganic phosphate, while a and 
(3 are the rates of the forward and backward reactions. The energy involved in the movement of 3 Na"*" and 2 ions 
against their electrochemical gradients is: 
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AwNa = -3e(t; - -WNa) (27) 
AwK = +2e{v - vk) , (28) 
where vk and ?;Na are given by equations (^ and (|^). This result is independant of the detailed interaction between 



ions, molecules and the ATPase enzyme. Therefore, the total work done in reaction ( Pq ) is: 

Aw — AwATP + AwNa + Awk 

= AwATP + e(3i;Na - 2wk - v) , (29) 

where Awatp is the energy associated with the breakdown of ATP. The ratio of a to is again given by equations 
( p^ ) and (p^). The energy available from ATP breakdown is much larger than the energy required to translocate 
the potassium and sodium ions at small negative membrane potentials. In practice, such a pump or motorized swing 
door will soon reach saturation. We therefore choose the rate of the forward reaction to be constant, resembling the 
maximum possible speed of the swing door in that direction, and obtain: 

a = A (30) 
" AwATP + e(3i;Na - 2wk - v) 



/3 = Aexp 

where A is a constant. The net pump current for a cell with M pumps can then be written: 



(31) 



a\ 1 /i I - 2wk + 3wNa + "atp I \ .„„^ 
«NaK = eM{a ~ P)^ fcNaK 1 " Gxp < ^ ^ , (32) 



where fcNaK = eAAf, and v^tv — Awatp/c 



D. Na+,Ca^+ Exchanger 



To maintain a steady state for the intracellular calcium concentration in for example heart cells, the amount of 
calcium that enters the cell via ionic channels must be extruded. The Na"*", Ca^"*" exchanger is the major mechanism 
responsible for achieving a balance between calcium entry and extrusion in oscillating cells. We assume that the rates 
of the exchange reaction (MuUins, 1977) given by: 



3Na^ 



Caf+ ^ 3Na+ + Ca2+ 
1^1 



(33) 



are governed largely by the electrochemical gradients for sodium and calcium, together with the membrane potential. 
In other words, the energy produced when 3 extracellular sodium ions take the elevator down into the cytosol is used 
to elevate one calcium ion up from the cytosol into the extracellular space: 



AlUNa = +3e(w - UNa) 

Awca = -2e(w - UCa) , 

where wca and WNa are given by equations (^ and (H). The total work done in reaction (^3|) is: 

Aw AwNa + Awca = e{v - 3uNa + 2uCa) • 



(34) 
(35) 

(36) 



The ratio of a to /3 in equation (p3| ) is given by equations ( [Tsl) and (|T^). Since Aw = e{v — 3uNa + 2wca) in a cell 
usually oscillates around zero, the most natural choice of a and /3 is: 



a = Aexp 
P = Aexp 



e{v - 3wNa + 2uca) 



2kT 

e{v - 3wNa + 2uca) 



2kT 



(37) 
(38) 



where again we make the ad hoc assumption that A is a constant. For a cell with N exchangers the net current is 
then: 



«NaCa = -eN (a - /?) = k 



NaCa 



sinh ■ 



V - 3DNa + 2t;ca 
2fcT 



(39) 



where /cNaCa = 2eXN. 
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E. Membrane Voltage 



We assume that the electrical activity of a cell is described by the five currents discussed above, and that all the 

other currents (Boyett, 1996) are of minor importance. The standard differential equation for the voltage, and the 
conservation laws for intracellular ionic concentrations are then: 

dv 1 . . . . 

= («K + «Ca + «Na + «NaCa + «NaKj , (40) 

IfKl- - ^'^"^ ~ f41) 

dt^ FV ' ^ ' 

J 2zNaCa^ (42) 

dV ^ 2FV ' ^ ' 

d I — JNa — S^NaK — SzNaCa 

^[Nah = — , (43) 

where C is cell capacitance, F is Faraday's constant and we assume that the cell volume V is constant. Solving 
equation ( pi] ) for ik, equation ( p2[ ) for ica, and equation ( p3[ ) for ZNa, we obtain: 

iK = -FV^\K]i + 2i^,K, (44) 

ZCa = -2FF^ [Ca]i + 2iNaCa , (45) 

dt 

«Na = -i^V^^iNa]! - 3iNaK - 3lNaCa • (46) 

Inserting this into equation ( ^0| ) we obtain: 



since all of the currents cancel. This equation can also be written as: 



Jt{^~^ {[K]i + 2[Ca]i + [Na]i}) = . (48) 

Integrating gives: 

FV 

— ([K]i + 2[Ca]i + [Na]i) = c , (49) 

where the integration constant c has to be determined. Given that the voltage across a capacitor is zero when the net 
charge difference is zero, choosing the integration constant as: 

FV 

c = -— {[K]c + 2[Ca]e + [Na] J , (50) 

gives: 

'' = ^ {(Mi - [K]c) + 2([Ca]i - [Ca]e) + ([Na]; - [Na]e)} , (51) 

so that the voltage is zero when there is no charge gradient. Equation (|5^) is simply the relation between the electric 
potential and charge of a capacitor, which is in fact the origin of equation (^). Thus it is completely general and 
independent of the number of membrane currents in a model. The meaning of equation ( |5l] ) is: 

The voltage across the membrane of a cell is caused by, and is directly proportional to, the surplus of charge 
inside the cell, 

but it is only valid if we assume that the intracellular concentration of anions is equal to the extracellular concentration 
of anions. Since the cations are balanced by the anions in any solution, the intracellular/extracellular concentration 
of anions must be close to the intracellular/extracellular concentration of cations. 
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Since equation (O) is the explicit integral of equation (^), it should be used instead of equation (40) (or the equiv- 
alent of equation (pO|)) in any model. The differential equation (^0|) is needed only in models where the intracellular 
ionic concentrations are not tracked individually (for example, the Hodgkin-Huxley equations (1952)). 

There is a significant difference between equations (^0|) and (jsj). Here is an analogy: There are two different ways to 
calculate how many ions are inside a cell. The first method counts every ion enetering or leaving (equation (^)), while 
the second method simultaneously counts all the ions inside the cell, perhaps by means of an ion detector (equation 
(pT|)). Both methods will correctly give the variation in the number of ions in the cell. However, the observer of ions 
entering and leaving will only observe variations in the number of ions. If he wants to know the actual number of 
ions in the cell, he must make an initial guess of the number of ions already inside. Because his guess may differ 
significantly from the actual ion number (as maesured by the ion detector) indicates, the results from the two methods 
may be contradictory. 



F. Osmotic Pressure 



In this section we demonstrate that variations in the osmotic pressure are a natural consequence of our model. In 
fact, van't Hoff's equation for the osmotic pressure can be derived by tracking the flow of energy under the assumption 
of constant temperature and volume. 

From elementary physics the work W done by a force J- to move an object from X = A to X = B is given by: 



W = / TdX ^ T{B -A)^ ^-r^^ = / ^^^^ ' 
Ja Jo dt Jq 



(52) 



where V = ^ is the velocity of the object, and X{0) — A, X{t) = B. When ions move across the membrane of a cell 
the transmembrane voltage v is like a force and the current i is like a velocity, thus: 

[ vidt. (53) 
Jo 

If we take into account the reversal potentials in our model this can be written: 

iK{v - Vk) + ica{v - VCa) + iNa{v - VNa) 

Jo 

+ iNaCa(w - 3wNa + 2uCa) + iNaK(w + 2wk " 3wNa)rft ■ (54) 



Using equations (40), ([44|), ( [45| ) and (^6|) the currents cancel and we obtain: 

fv .[K]i .[Ca], .[Na]i 

W^-C vdv + FV VKd{[K]i) + 2FV vcadi[Ga]i) + FV W([Na]i) (55) 

"'O ^[K]<3 "'[Ca]c "'[Na]o 

When solving these integrals it is important to remember that the reversal potentials vk, ^ca and WNa given by 
equations (H) and (^, are dependent on the integration variables [K]i, [Ca]i and [Na]i. Using: 



\n{(p)d(p ^ (f)ln{(p) - (f> , (56) 



the integral becomes: 



W = -Icv' - RTv([K], In(j^) + [Ca]i In(i^) + [Na]i ln( 



2 [K]c' ' ' MCa]e' ' ' MNa]e'J (57) 

- RTV{[K], - [K]i + [Na]e - [Na]i + [Ca]c - [Ca];} . 

This is the total work associated with the currents. The pump work and work due to the exchanger and channel 
currents have opposite signs since the pump moves ions against their electrochemical gradients, while energy is lost 
by the downhill ionic fluxes through the exchanger and the channels. Since energy cannot be created or destroyed, 
we define the potential energy P of the cell so that the total energy E is zero (this is the most common definition of 
potential energy): 
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E = P + W = , (58) 

thus, 

P ^ lev + RTv{m. H&) + [C], in(gli) + INal, j^J)} 

+ i?TF {[K]e - [K]i + [Na]e - [Na]i + [Ca]e - [Ca];} . 

The first term is the electrical potential energy of a capacitor, the second term is the potential energy associated with 
the concentration gradients, and the third term is the potential energy associated with an osmotic pressure gradient 
(caused by the concentration gradients). Rewriting equation (p8[): 

t 

iNaK(w + 2UK - 3uNa)rfi = (60) 

Icv' (61) 
RTv[m ln(|^) + [Ca]. ln(a) + [Na]. ln(»)} (62) 
RTV {[K]e - [K]i + [Na]e - [Na]i + [Ca]e - [Ca];} (63) 

iK{v - Vk) + «Ca(w - WCa) + iNa(w " ^Na) + iNaCa(w " 3wNa + 2vcii)dt , (64) 

we see that the pump (equation (|60|)) produces: 

1. a transmembrane voltage gradient (equation (^l])), 

2. a transmembrane concentration gradient (equation (|62[)), 

3. a transmembrane pressure gradient (equation (pq)), and 



4. downhill fluxes through the exchanger and the channels (equation (64)). 

In 1886 van't HofF noticed that the behaviour of solutes in dilute solutions resembles the behaviour of a perfect gas: 

The pressure which a gas exerts when a given number of molecules are distributed in a given volume is equally great as 
the osmotic pressure, which under the same conditions would be produced my most solutes when they are dissolved in an 
arbitrary solvent. 

Therefore we can use the ideal gas law p — [S]RT to calulate the osmotic pressure across the cell membrane: 

^ - i?T { [K]i - [K]e + [Na]i - [Na]e + [Ca]i - [Ca]c} . (65) 

This is the van't Hoff equation for the osmotic pressure (tt) across a solute impermeable barrier separating two ideal 
dilute solutions. This is nothing but equation ( |63|) divided by the cell volume V. 

Since we did not have the osmotic pressure in mind when we made the theory, our theory is a perfect example of 
one of Richard P. Feynman's astute observations: 

When you have put a lot of ideas together to make an elaborate theory, you want to make sure, when explaining what it 
fits, that those things it fits are not just the things that gave you the idea of the theory; but that the finished theory makes 
something else come out right, in addition. 

Our theory indeed has "something else [that] comes out right" - the osmotic pressure variations. 
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III. MODEL PARAMETERS 



A mathematical model of the membrane potential has been derived; equations 
equilibrium potentials, equations ( p2|), ( |2^ ) and (^5|) the ionic currents, equations (fe 
the pump currents, equations (41), (|42|) and ( ^ ) the ionic concentrations, equation 

finally, equation (^5|) the osmotic pressure across the cell membrane. The model has 6 dimensions with the variables 
X, /, h, [K]i, [Ca]i and [Na];. 



, (H) and represent the 
and the exchanger and 
the membrane voltage, and 



2[K]. 



3[Na] 




I 



[Ca ] 

3[Na] 3[Na] 



ATP 



ADP 




FIG. 1. Schematic diagram of the simplified cell model. The three ionic currents of potassium, calcium and sodium, the 
exchanger and the pump. The "moveable" gating doors of the ionic channels illustrate that the concentration gradient of 
potassium is the opposite to the concentration gradients of calcium and sodium. 

Fig. 1 shows a schematic diagram of the model, with the ionic channels of potassium, calcium and sodium, the 
Na^,Ca'^^ exchanger and the Na''',K+ pump. The "moveable" gating doors of the ionic channels illustrate that the 
concentration gradient of potassium is the opposite of the concentration gradients of calcium and sodium. 

On the basis of distinct biophysical and pharmacological properties, cardiac calcium currents have been classified 
into a long lasting type (ica.L) and a transient type (ica.x), while cardiac delayed rectifier potassium currents have 
been classified into a rapid type (iKr) and a slow type (zks)- In our model we assume that ica.T and ixs are of minor 
importance; i.e. when we talk about ica we mean ica.L, and when we talk about «k we mean zki-- 

First, we want to justify the presence of the term kT/2e in equations (|2^), (24) and ([25|). This corre- 
sponds to a slope factor for the activation and inactivation curves of kT/Ae ~ 6.68 mV at 37°C. The ob- 
served slope factors are 7.4 mV for activation of zk (Shibasaki, 1987), 6.6 mV for activation of ica (Hagi- 
wara et ai, 1988), 6.0 mV for inactivation of ica (Hagiwara et ai, 1988), 6.0 mV for activation of ZNa (Mu- 
ramatsu et ai, 1996), and, finally, 6.4 mV for inactivation of iNa (Muramatsu et ai, 1996). Hence, we see 
that kT/Ae^ corresponding to a gating charge of g « 4e, is an excellent approximation. Now, we would 
like to distinguish between the fundamental physical constants (table ||), the experimentally observed con- 
stants (table the adjustable parameters (table [ll^ and the initial conditions (table ^V[ ) in the model. 



TABLE L 



Fundamental Physical Constants 



Name 


Value 


Unit 


k 


1.38065812 • 10"^" 


mJ/K 


e 


1.6021773349 • 10"^^ 


C 


F 


96485.30929 


C/mol 


T 


310.15 


K 


R = kF/e 


8314.511935 


J/kmol K 



The half-activation and inactivation potentials in the model (wx, Wd, wf, and Vh) are based on the ex- 
periments of Shibasaki (1987), Hagiwara et al. (1988) and Muramatsu et al. (1996), and we use a value of 
fATP that gives a reversal potential for the sodium pump in good agreement with the experiments of Sakai et 
al. (1996). The maximum time constants in these experiments were 203ms for activation of (Shibasaki, 
1987), 225 ms for inactivation of ica (Hagiwara et al., 1988) and 174ms for inactivation of iNa (Muramatsu et 
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UjL. ^ A.U tJyJ I • 


In the model, however, we combine these and 


use a maximum time constant of 200 ms for 


both Tk, 


TCa and TNa 


. Finally, we use typical values for cell volume, cell capacitance, and extracellular ionic concentrations. 


TA RT TT' TT 




Observed Constants 




Name 


Value 




Unit 


[K]c 






mi\/r 

niivi 


[Ca]e 


9 




mivj. 


[Na]e 


1 AC\ 




llliVi 


V 


xu 




1 iim^ 
lU fJ/lLL 


C 


4:1 




pi- 


Wx 


-25.1 




mV 


Vd 


-6.6 




mV 


Vi 


-25.0 




mV 


Vm 


-41.4 




mV 


Vh 


-91.0 




mV 


vatp 


-450 




mV 


T 


200 




ms 



The density of ionic channels, exchangers and pumps (i.e. gca, ffNa, Sk, ^NaK and fcNaCa) can vary signifi- 
cantly from cell to cell. In order to reproduce the action potentials recorded in a spontaneously beating adult 
rabbit sinoatrial node cell with normal Tyrode solution as external super fusing solution (Fig. 7 A. in Barus- 



cotti et al. (1996)), we fit the adjustable parameters (table III) and the initial conditions (table |IV| ) numeri- 
cally. More details are given in (Endresen, 1997a). Many combinations of gca, gNa, ffK, ^NaK and /cNaCa re- 
sulted in good approximations to the experimentally recorded waveform, from which we conclude that different 
cells can produce the same action potential although they have a different mixture of ionic channels, exchang- 



ers and pumps. The final set of parameters presented in table [II are based on the choice /sNaK = 12.2 pA 



TABLE III. Adjustable Parameters 

Name Value Unit 

g^^ 9.29045 nS~ 

SNa 253.94203 nS 

QK 0.70302 nS 

fcNaCa 8181.31568 pA 

fcNaK 12.20000 pA 



TABLE IV. Initial Conditions 



Name Value Unit 

xo - 

/o 1 

ho - 

[K]ig 130.880955 mM 

[Ca]ig 0.000790 mM 

[Na]io 18.514880 mM 
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IV. SIMULATION RESULTS 



The six differential equations in the model were solved numerically using a fifth-order Runge-Kutta method with 
variable steplength, more details are given in (Endresen, 1997b). Fig. 2 shows the five membrane currents and Fig. 
3 shows the recorded and simulated action potentiald together with the intracellular ionic concentrations. These 
computations used the initial conditions stated in table IV. 

Cells must generate their membrane potential by actively transporting ions against their concentration gradients. 
To examine this process in our model, we ran a simulation starting with equal intracellular and extracellular ionic 
concentrations: 



[K]i = 5.4 mM , 

[Ca]i = 2mM, (66) 
[Na]i = 140 niM . 



400 



300 - 
200 - 




-300 - 



-400 



-500 







Us) 



FIG. 2. Membrane currents in a simple ionic model of a rabbit sinoatrial node cell. We show the outward delayed rectifying 
potassium current (ik), the inward calcium current (ica), the inward sodium current (iNa), the sodium calcium exchange 
current (iNaCa) and the sodium potassium pump current (iNaK). These computations used the initial conditions in table IV. 



The results are presented in Fig. 4. After approximately 1500 seconds (25 minutes) the system appears to reach a 
stable fixed point with the following intracellular ionic concentrations: 



[K]i ^ 115.842881 mM , 
[Ca]i = 4.485016- 10"^ mM 
[Na]i = 33.548671 mM . 



(67) 



At 2500 seconds we kick the cell with a 20 pA pulse of potassium ions with duration 50 ms to see if it wants to oscillate. 
It surely does, and the oscillations soon become almost identical to the original oscillations present in Fig. 2 and 3 
(you can not see this from Fig. 4 since the time scale is much to big). The values for the concentrations only differ 
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slightly from the initial conditions in table IV , as indicated by the final concentrations calculated by the program (at 
5000 seconds): 



[K]i = 131.075490 mM , 

[Ca]i = 6.827191 • 10"^ mM , (68) 
[Na]i = 18.320693 mM . 




FIG. 3. Membrane potential, intracellular ionic concentrations and osmotic pressure of a rabbit sinoatrial node cell, (a) 
Model-generated (solid line) and experimentally recorded (dashed line) action potential waveform, (b) potassium concentration 
[K]i, (c) calcium concentration [Ca]i, (d) sodium concentration [Na]i and (e) the osmotic pressure tt across the cell membrane. 
These computations used the initial conditions in table [iv] . 



The trajectory reaches the limit cycle at 2560 seconds, and there is no drift in the intracellular ionic concentrations 
from this point to the end of the simulation (at 5000 seconds). The long time simulation in Fig. 4 is a numerical 
indication that the oscillation in Fig. 2 and 3 is indeed a stable limit cycle. However, the model also has a stable 
fixed point given by (|67|). 



V. DISCUSSION 

We have presented a simple model for the cells of the rabbit sinoatrial node. Our model involves only Na''',K+, 
and Ca^'^ ions, their respective channels, the Na^,Ca^^ exchanger, and the Na^,K+ pump. The equations were 
derived using basic physical principles and conservation laws. Since the only source of energy in our model is the 
sodium pump, we could easily track the flow of energy in our system. We showed that the pump works to generate 
a transmembrane voltage, pressure gradient and ionic concentration gradients. Our equations also accounted for 
the energy lost due to downhill ionic fluxes through the exchanger and channels. A prediction of osmotic pressure 
variations was a novel result of our energy analysis. 
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The intracellular ionic concentrations are dynamic variables in our model, governed by the conservation equations 
(p]), p^), and ( ^). This allowed us to replace the standard differential equation for the voltage ( ^0| ) with the 
algebraic equation (|5l|) . Although a number of other ionic models also keep track of intracellular ionic concentrations 
(see Wilders (1993)), we are unaware of any other model using an algebraic equation for the membrane potential. 
Models that use the standard voltage differential equation ( ^ ) have a superfluous extra dimension. Furthermore, 
the initial conditions for this extra differential equation cannot be chosen independently of the initial conditions of 
the conservation equations (p|), (^2|), and (^) - otherwise, the computed membrane potential will be erroneous. 
For these reasons, we suggest that our algebraic expression for the membrane potential should replace the standard 
voltage differential equation in models where intracellular ionic concentrations are dynamic variables. 
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FIG. 4. Long time simulation showing the membrane potential, intracellular ionic concentrations and osmotic pressure 
starting with equal intracellular and extracellular concentrations: [K]i = 5.4 mM, [Ca]i = 2mM and [Na]i = 140 mM. The cell 
is given a 20 pA pulse of potassium with duration 50 ms after 2500 seconds of the simulation, (a) Membrane potential v, (b) 
potassium concentration [K]i, (c) calcium concentration [Ca]i, (d) sodium concentration [Na]i and (e) the osmotic pressure tt 
across the cell membrane. 



Our model does not include funny current (if), ATP sensitive channels, stretch-activated channels, or other ion 
channels that may be important (Boyett, 1996). We also ignored the effect of calcium uptake and release from 
the sarcoplasmatic reticulum, and the effect of anions like chloride. We have assumed that the ionic channels are 
governed by a Markov process. However, since transmembrane proteins can have a large number of conformational 
states (Liebovitch, 1995), perhaps a fractal model for the gating would be a better approach. We assumed that 
the maximum of the activation/inactivation time constant occurs at the same voltage as the inflection point of the 
sigmoidal steady state activation/inactivation curve. Also, we have assumed that the cell volume is constant. While 
such assumptions reduce the number of parameters in the model, they may also result in discrepancies with experiment. 
A natural extension of our model would include a variable cell volume. 

Finally, we would like to point out that our model is based on experiments where some were conducted at room 
temperature (22-24°C) (Baruscotti et ai, 1996; Muramatsu et ai, 1996), while others were performed at 37°C 
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(Shibasaki, 1987; Hagiwara et ai, 1988; Sakai et al, 1996). It is not clear what affect this inconsistency had on the 
accuracy of our model. 

The values of the parameters gca, 5Na, 5k, ^NaK and /cNaCa, given in table [II, are only an estimate of the actual 
physiological parameters. We did not systematically study the dynamics of the model for different parameters. How- 
ever, interested readers can download the source code for the model (in both UNIX and Windows NT environments) 
from the following internet site: 



http://www.physio.mcgill.ca/guevaralab/singlecell.htm 



We hope that future experiments will help us to discriminate between the different parameter sets that reproduce the 
experimentally recorded action potentials. 
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